/**
 * @integralImageConvolution.cpp
 * @brief Questo header contiene il corpo delle funzioni per l'applicazione
 * di un filtro gaussiano sfruttando la tecnica dell'integral image
 *
 * 
 *
 * @date 22/11/2012
 */
#include "integralImageConvolution.h"

void calcIntegralImage(const Mat_<uchar>& image, Mat_<int>& intImg){

    intImg = Mat_<int>::zeros(image.rows, image.cols);

    // Sum along rows
    for(int r=0; r < intImg.rows; ++r){
        // Initialize the first column of integral image
        intImg(r, 0) = image(r, 0);

        for(int c=0; c < intImg.cols-1; ++c)
            intImg(r, c+1) = intImg(r,c) + image(r, c+1);
    }

    // Sum along columns
    for(int c=0; c < intImg.cols; ++c)
        for(int r=0; r < intImg.rows-1; ++r)
            intImg(r+1, c) = intImg(r,c) + intImg(r+1, c);

}

void medianFilterIntegImage(const Mat_<int>& intImg, Mat_<uchar>& imgOut, const int kern_dim){

    // kernel must be at least 3x3
    if(kern_dim < 3) return;

    // Calculate half dim of kernel matrix (integer division)
    int k_w = kern_dim/2;

    int n_kern_elem = kern_dim*kern_dim;

    for(int r=k_w+1; r<intImg.rows-k_w; ++r){
        for(int c=k_w+1; c<intImg.cols-k_w; ++c){

            //Calc area of kernel window centered in current pixel
            int area = intImg(r+k_w, c+k_w) - intImg(r+k_w, c-k_w-1) - intImg(r-k_w-1, c+k_w) + intImg(r-k_w-1, c-k_w-1);

            //Assign median value to output image
            imgOut(r,c) = area/n_kern_elem;
        }
    }
}

void integralImgGaussianConv(const Mat_<uchar> &imgIn, Mat_<uchar>& imgOut, const double sigma){

    //imgOut = Mat::zeros(imgIn.size(), CV_8U);

    double TwelveSigmaQuad = 12*sigma*sigma;
    // Ideal kernel width for building gaussian filter with var==sigma
    int Wid = floor(sqrt((TwelveSigmaQuad/NPASSES) + 1));

    // Calculate width of the two kernels to apply
    int Wl;
    (Wid%2) > 0 ? Wl = Wid : Wl = Wid-1;
    int Wu = Wl + 2;

    //Find num times we should apply Wl kernel
    int m = round( (TwelveSigmaQuad - NPASSES*Wl*Wl - 4*NPASSES*Wl - 3*NPASSES)/(-4*Wl - 4) );

    imgIn.copyTo(imgOut);

    Mat_<int> intImg;

    if(Wl >= 3){    //If Wl kernel dim is at least 3
        for(int i=0; i < m; ++i){
            calcIntegralImage(imgOut, intImg);
            medianFilterIntegImage(intImg, imgOut, Wl);
        }
    }

    for(int i=0; i < NPASSES-m; ++i){
        calcIntegralImage(imgOut, intImg);
        medianFilterIntegImage(intImg, imgOut, Wu);
    }
    intImg.release();
}

// ----------------------- IMPLEMENTATION WITH POINTERS --------------------------------

void calcIntegralImagePoint(const Mat_<uchar>& image, Mat_<int>& intImg){

    int img_rows = image.rows;
    int img_cols = image.cols;

    intImg = Mat_<int>::zeros(img_rows, img_cols);

    //Take pointers to Mat's data
    uchar* img = image.data;
    int* intimg = (int*)intImg.data;

    // Sum along rows
    for(int r=0; r < img_rows; ++r){

        int row_first_px = r*img_cols;
        // Initialize the first column of integral image
        *(intimg + r*img_cols) = *(img + r*img_cols);

        for(int c=0; c < img_cols-1; ++c)
            *(intimg + row_first_px + c+1) = *(intimg + row_first_px +c) + *(img + row_first_px + c+1);

    }

    // Sum along columns
    for(int c=0; c < intImg.cols; ++c)
        for(int r=0; r < intImg.rows-1; ++r){

            int row_first_px = r*img_cols;
            int rowplus1 = (r+1)*img_cols;

            *(intimg + rowplus1 +c) = *(intimg + row_first_px +c) + *(intimg + rowplus1 +c);
        }
}

void medianFilterIntegImagePoint(const Mat_<int>& intImg, Mat_<uchar>& imgOut, const int kern_dim){

    // kernel must be at least 3x3
    if(kern_dim < 3) return;

    int* intimg = (int*) intImg.data;
    uchar* imgout = imgOut.data;

    int img_rows = intImg.rows;
    int img_cols = intImg.cols;

    // Calculate half dim of kernel matrix (integer division)
    int k_w = kern_dim/2;

    int n_kern_elem = kern_dim*kern_dim;

    for(int r=k_w+1; r<img_rows-k_w; ++r){
        for(int c=k_w+1; c<img_cols-k_w; ++c){

            int row_first_px = r*img_cols;
            int rowpluskw = (r+k_w)*img_cols;
            int row_minkw_min1 = (r-k_w-1)*img_cols;

            //Calc area of kernel window centered in current pixel
            int area = *(intimg + rowpluskw + c+k_w) - *(intimg + rowpluskw + c-k_w-1) -
                    *(intimg + row_minkw_min1 +c+k_w) + *(intimg + row_minkw_min1 + c-k_w-1);

            //Assign median value to output image
            *(imgout + row_first_px + c) = area/n_kern_elem;
        }
    }
}

void integralImgGaussianConvPoint(const Mat_<uchar> &imgIn, Mat_<uchar>& imgOut, const double sigma){

    //imgOut = Mat::zeros(imgIn.size(), CV_8U);

    double TwelveSigmaQuad = 12*sigma*sigma;
    // Ideal kernel width for building gaussian filter with var==sigma
    int Wid = floor(sqrt((TwelveSigmaQuad/NPASSES) + 1));

    // Calculate width of the two kernels to apply
    int Wl;
    (Wid%2) > 0 ? Wl = Wid : Wl = Wid-1;
    int Wu = Wl + 2;

    //Find num times we should apply Wl kernel
    int m = round( (TwelveSigmaQuad - NPASSES*Wl*Wl - 4*NPASSES*Wl - 3*NPASSES)/(-4*Wl - 4) );

    imgIn.copyTo(imgOut);

    Mat_<int> intImg;

    if(Wl >= 3){    //If Wl kernel dim is at least 3
        for(int i=0; i < m; ++i){
            calcIntegralImagePoint(imgOut, intImg);
            medianFilterIntegImagePoint(intImg, imgOut, Wl);
        }
    }

    for(int i=0; i < NPASSES-m; ++i){
        calcIntegralImagePoint(imgOut, intImg);
        medianFilterIntegImagePoint(intImg, imgOut, Wu);
    }
    intImg.release();
}

double round(double x){

	return floor(x+0.5);
}